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The  quasi-statie  mieromeehanieal  damage  mechanics  originally  formulated  by  Ashby  and  Sammis  (PAGEOPH, 
1990)  has  been  expanded  in  three  important  ways:  1)  An  energy  density  function  has  been  derived  that  allows  a 
self-consistent  inclusion  of  the  effects  of  dynamic  damage  evolution  on  the  elastic  and  anelastie  response,  2) 
whereas  the  Ashby/Sammis  model  was  only  applicable  to  the  ease  where  the  initial  cracks  are  all  parallel  and  the 
same  size,  we  can  now  include  a  specified  distribution  of  initial  eraek  sizes  with  random  azimuthal  orientation  about 
the  loading  axis,  and  3)  we  allow  for  yielding  of  the  weaker  minerals  in  granite  at  relatively  low  stress  levels  which 
produce  strong  nonlinearity  in  the  failure  envelope  and  stress-strain  eurve.  This  new  damage  meehames  is  validated 
using  published  experimental  data  for  Westerly  granite.  In  collaboration  with  the  mechanics  group  at  U.C.  Santa 
Barbara,  we  have  built  our  new  damage  mechanics  into  the  commercial  finite  element  program  ABAQUS  and  used 
it  to  simulate  an  explosion  in  damaged  rock.  We  find  that  an  explosion  in  a  pre-stressed  medium  or  one  that  contains 
anisotropic  initial  damage  generates  strong  5-wave  radiation  in  agreement  with  previous  simpler  calculations  by 
Johnson  and  Sammis  (PAGEOPH,  2001)  and  experimental  observations  in  damaged  photoelastie  plates. 
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OBJECTIVES 

The  overall  objective  of  this  research  program  has  been  to  understand  and  quantify  the  extent  of  fracture  damage  in 
the  non-linear  source  region  of  an  underground  nuclear  explosion  and  to  assess  its  effect  on  the  radiation  of  seismic 
energy  to  the  far  field.  Of  particular  interest  is  the  generation  of  high-frequency  shear  waves  that  may  affect  regional 
discrimination  and  yield  estimation.  The  specific  objective  of  this  year’s  research  has  been  to  build  on  the  work  of 
Deshpande  and  Evans  (2008),  to  extend  the  micromechanical  damage  mechanics  originally  formulated  by  Ashby 
and  Sammis  (1990),  and  to  incorporate  it  in  numerical  models  that  simulate  underground  explosions. 


RESEARCH  ACCOMPLISHED 

Our  major  accomplishment  this  year  has  been  to  make  three  important  improvements  in  the  micromechanical 
damage  mechanics  developed  by  Ashby  and  Sammis  (1990)  and  expanded  by  Deshpande  and  Evans  (2008).  First, 
we  dropped  the  assumption,  implicit  in  both  formulations,  that  the  growing  wing  cracks  are  all  parallel  and  all  the 
same  size,  which  strongly  affects  the  interaction  between  the  growing  wing  cracks  that  leads  to  failure.  Second,  we 
calculated  the  change  in  elastic  energy  density  associated  with  the  growing  wing  cracks  in  a  less  approximate  way 
than  that  used  by  Deshpande  and  Evans  (2008),  which  has  led  to  a  more  accurate  stress-strain  simulation.  Third,  we 
account  for  the  fact  that  granite  is  multimineralic  and  that  yielding  in  the  weaker  minerals  begins  at  relatively  low 
stress  levels. 

We  have  built  this  new  mechanics  into  the  dynamic  finite  element  code  ABAQUS  and  produced  some  initial 
simulations  of  damage  evolution  during  an  underground  explosion.  Finally,  we  have  begun  experiments  in  which  a 
hyper-velocity  impact  is  used  to  simulate  an  explosion.  A  hyper-velocity  impact  creates  a  much  higher  energy 
density  at  the  source  than  does  an  exploding  wire,  which  allows  us  to  better  scale  the  damage  process  from  the  field 
to  the  lab. 

Improvmcnts  to  the  Ashby  Sammis  (1900)  Micromechanical  Damage  Mechanics 

Ashby  and  Sammis  (1990)  (hereafter  referred  to  as  A&S)  formulated  a  micromechanical  damage  mechanics  based 
on  the  nucleation,  growth  and  interaction  of  tensile  “wing  cracks”  nucleated  at  the  tips  of  an  initial  distribution  of 
penny-shaped  microcracks.  The  basic  elements  of  their  model  are  illustrated  in  Figure  1. 

A&S  demonstrated  that  their  model  gives  an  adequate  description  of  the  failure  envelope  (oj  vs  cr3  at  failure)  at  low 
confining  pressures,  ay  for  a  w  ide  range  of  rocks  loaded  in  triaxial  compression  (o^  <  as  =  ay  where  compression  is 
taken  as  negative). 
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Figure  1.  Geometry  of  damage  mechanics 
model  (from  Deshpande  and  Evans, 
2008).  And  initial  distribution  of  flaws 
of  radius  a  are  loaded  in  compression. 
Sliding  on  the  initial  flaws  nucleates 
tensile  “wing  cracks”  at  their  tips  of 
length  /.  The  wedging  force  produced 
by  the  angle  cracks  creates  tension  on 
the  remaining  ligament  creating  a 
positive  feedback  that  leads  to 
instability  and  failure.  The  initial 
damage  is  defined  as  D0=^7r(aaf  Ny 

where  ;VV is  the  volume  density  of  initial 
flaws.  Damage  increases  as  the  wing 
cracks  grow  according  to 

D  =  ^  x(l  +  aaf  Nv  • 
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However,  the  failure  envelope  predicted  by  this  model  is  nearly  linear  while  the  observed  failure  envelope  for 
Westerly  granite  is  very  non-linear,  especially  at  high  eonfining  pressures  (see  Figures  2  and  3).  A&S  hypothesized 
that  this  nonlinearity  is  caused  by  a  transition  to  plastie  deformation  at  very  high  eonfining  pressures.  We  shall  show 
in  a  subsequent  seetion  that  the  broad  transition  zone  in  Fig.  3  is  due  to  the  multimineralie  composition  of  Westerly 
granite  where  the  weaker  minerals  yield  at  relatively  low  stress  levels. 


Figure  2.  Failure  envelope  for  Westerly 
granite  under  triaxial  loading 
<j\  <  <r?  =  cry  Sources  of  data  are 
labeled.  The  heavy  dashed  line 
shows  the  yield  stress  for  quartz. 
The  light  dashed  line  shows  the 
stress  state  at  which  the  initial  flaws 
nucleate  wing  cracks.  Note  that  no 
additional  damage  can  occur  for 
eonfining  stress  greater  than  about 
1000  MPa  (confining  pressure  is 
taken  as  negative). 
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Figure  3.  Enlarged  view  of  failure  envelope  in 
Fig.  2  near  <r3  =0.  The  dashed  curve  is 
the  theoretical  failure  envelope  when  all 
flaws  are  the  same  size.  The  variable  pis 
a  correction  for  a  distribution  or  crack 
orientations  defined  in  the  text.  Note 
that  the  departure  from  linearity  begins 
at  <jy  %  100  MPa.  The  key  question  is 
whether  this  deviation  from  linearity  is 
due  to  a  distribution  of  sizes  of  starter 
flaws  or  the  mineralic  heterogeneity  of 
Westerly  granite  where  the  weaker 
grains  begin  to  yield  at  significantly 


Johnson  and  Sammis  (2001)  used  the  A&S  damage  mechanics  to  model  the  1  ktehemieal  explosion  detonated  in 
September  1993  as  part  of  the  Non-Proliferation  Experiment  (NPE).  They  found  that  the  explosion  indueed  fraeture 
damage  to  a  distance  up  to  ten  times  the  eavity  radius.  More  importantly,  they  found  that  the  damage  itself  generated 
significant  “secondary  radiation”  of  P  and  S  waves  in  the  seismic  band.  However,  the  A&S  formulation  they  used  is 
incomplete.  Although  it  allows  an  estimate  of  the  evolving  damage  and  predicts  failure,  it  does  not  provide  an 
estimate  of  the  stress-strain  behavior  associated  w  ith  the  increasing  damage.  For  their  numerical  simulations, 
Johnson  and  Sammis  (2001)  used  an  empirieal  stress-strain  relation  borrowed  from  soil  meehanies,  whieh  was  not 
connected  to  the  damage  model. 

Deshpande  and  Evans  (2006)  (hereafter  referred  to  as  D&E)  extended  the  A&S  damage  meehanies  to  include 
stress-strain  behavior  by  deriving  the  internal  energy  density  as  a  function  of  stress,  whieh  they  then  differentiated  to 
find  the  strain.  They  wrote  the  energy  density  W as 


W  =  Wa  +  NyAiVl  +  NyMV2 


(1) 


559 


2009  Monitoring  Research  Review:  Ground-Based  Nuclear  Explosion  Monitoring  Technologies 


where  W a  is  the  elastic  energy  density,  a  IF,  is  the  change  in  energy  density  per  crack  due  to  sliding  on  the 
preexisting  cracks,  and  AW2  is  the  change  in  energy  density  per  crack  due  to  the  growth  of  wing  cracks.  The  axial 
strain  for  any  stress  state  (cr,,  o-3)  may  be  calculated  from  the  energy  density  as 

‘Wfa'T,)  (2) 

rb*i 

The  first  two  terms  in  the  energy  density  are  quadratic  in  cj\  and  thus  give  a  linear  stress-strain  curve  on 
differentiation.  The  second  term  reduces  the  effective  Young’s  modulus,  as  in  O’Connell  and  Budiansky  (1974). 
The  third  term  is  not  quadratic  and  produces  curvature  in  the  stress-strain  curve,  in  the  sense  of  the  observations  in 
Figures  2  and  3.  The  question  is,  how  much  curvature? 

The  perturbation  in  internal  energy  due  to  the  grow  th  of  wing  cracks  can  be  written 

AW2  =.—  Jujrdr  (3) 

E  n  __ 


where  K}  is  the  mode  I  stress  intensity  factor  at  the  tips  of  the  wing  cracks  and  /,  a ,  and  a  are  defined  in  Figure  1 . 
A&S  derive  an  analytical  expression  for  K{  that  can  be  integrated  to  find  AWV 

In  their  formulation,  D&E  simplify  the  integration  in  eqn.  (3)  by  assuming  the  wing  cracks  grow  at  a  fixed  value  of 
l/a  and  WTite 

AIV2=^  \ Kjrdr  (4) 

bo  o 

Since  K\  is  a  constant  for  fixed  l/a,  they  can  remove  it  from  the  integral  giving  a  very  simple  result.  However,  we 
have  shown  (Sammis  et  al.,  2009)  that  this  simplification  is  too  extreme  in  that  AW2  becomes  independent  of  the 
evolving  damage.  Hence  in  the  D&E  approximation,  wing-crack  growth  produces  zero  additional  strain  when  the 
strain  is  calculated  as  =cW  Ida ,  according  to  equation  (2). 


However,  Sammis  et  al.  (2009)  point  out  that  Kf  =Kk  (its  critical  value)  during  quasi-static  crack  growth  so  it  can  be 
removed  from  the  integral  in  eqn.  (3)  giving  an  equally  simple  result 


AW  2  =■  — 


3  KlD„ 


4  Enaa 


(  D  ^  3 


\E>0J 


-1 


4  3 

where  D()  =  —  x(aa)  Nv  is  the  initial  damage  associated  with  the  angled  cracks, 

damage,  and  Ny  is  the  volume  density  of  initial  flaws.  Since  D  increase  with  cr, 
contributed  by  the  growing  wing  cracks  is  non-zero. 


(5) 

4  3  . 

D()  -  —  n(l+  aa)  Ny  is  the  current 
up  to  failure,  the  nonlinear  strain 


Even  with  this  correction,  we  found  that  the  damage  model  fails  to  predict  the  observed  curvature  in  the  failure 
envelope.  Figure  3  compares  the  theoretical  and  observed  failure  envelopes  which  diverge  at  confining  pressures 
above  about  100  MPa.  Figure  4  shows  the  nonlinear  portion  of  the  stress  strain  curve  for  Westerly  granite  at  a 
confining  pressure  of  200  MPa  from  Lockner  (1998).  The  extra  strain  associated  with  the  growth  of  w  ing  cracks  (the 
only  nonlinear  contribution  in  the  model)  is  more  than  an  order  of  magnitude  smaller  than  the  measured 
nonlinearity. 
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Figure  4.  Non-linear  strain  in 

Westerly  granite  at  a  confining 
pressure  of  200  MPa  (redrawn 
from  Lockner,  1998).  Note  that 
significant  non-linear  strain 
begins  at  an  axial  stress  equal 
to  about  half  of  the  failure 
stress  and  that  the  non-linear 
portion  of  the  strain  at  failure 
is  about  0.003.  The  non-linear 
strain  associated  with  the 
growth  of  wing  cracks  (solid 
circles)  is  about  one  order  of 
magnitude  smaller. 
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Figure  5.  Same  as  Fig.4  but  with 
logarithmic  nonlinear 
strain  and  differential 
stress.  The  straight  line 
implies  that  the  stress- 
strain  curve  is  of  the  form 

a\  ~<7y- A  log£w,„/<wtir  +  B 


0.0001  0.001  0.01 

Nbnlinear  Strain 

By  dropping  two  simplifying  assumptions  made  by  both  A&S  and  D&E  the  damage  mechanics  can  produce  more 
nonlinearity  in  the  failure  envelope:  1)  we  no  longer  assume  that  all  the  initial  flaws  are  parallel  and  2)  we  no  longer 
assume  that  all  flaws  are  the  same  size.  We  can  drop  the  assumption  that  all  flaws  are  parallel  by  modifying  the 
interaction  term.  We  still  assume  that  all  flaws  are  all  at  an  angle  a  to  the  X/  loading  axis  but  allow  the  Xi  and  xj 
components  of  their  normal  vectors  to  be  randomly  distributed.  The  effect  is  to  reduce  both  the  wedging  force  and 
average  eross-seetion  area  supporting  the  interaction.  We  drop  the  seeond  assumption  that  all  flaws  are  the  same  size 
by  simply  keeping  track  of  the  growth  and  interaction  of  a  distribution  of  initial  flaw  sizes.  We  ehose  the  distribution 
measured  in  Westerly  granite  by  Hadley  (1976)  shown  in  Figure  6.  This  is  a  bimodal  distribution  having  a 
population  of  flaws  at  the  grain  boundary  scale  (about  1  mm)  and  a  seeond  power-law  (fractal)  distribution  of 
microcracks  at  the  10-100  micron  scale. 
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Figure  6.  Initial  flaw  distribution  in  Westerly 
granite  measured  by  Hadley  (1976). 
This  is  a  bimodal  distribution  with 
flaws  of  radius  a  near  the  size  of  the 
grain  boundaries  (about  0.5  mm) 
and  a  second  fractal  population  of 
microcracks  with  radii  less  than  100 
microns.  The  point  labeled 
“Bimodal  Density”  is  the  density  of 
smaller  flaws  required  to  produce 
the  observed  curvature  in  the 
failure  envelope.  However,  even 
with  this  larger  density,  the  stress- 
strain  curve  does  not  have  enough 
nonlinearity. 


Although  both  corrections  produced  more  damage  prior  to  failure,  the  total  strain  at  failure  (as  shown  by  the  closed 
circles  in  Fig.  4)  was  still  only  about  5xl0‘4,  a  full  order  of  magnitude  smaller  than  the  measured  value.  These 
corrections  also  failed  to  produce  the  observed  curvature  in  the  failure  envelope  in  Figures  2  and  3.  Even  with  the 
extra  damage,  the  failure  envelope  was  nearly  linear.  We  could  simulate  the  observed  curvature  in  the  failure 
envelope  by  increasing  the  density  of  microfractures  relative  to  the  grain  boundaries  (closed  circle  labeled  “bimodal 
density'’  in  Fig.  6),  but  this  required  micro-fracture  densities  two  orders  of  magnitude  larger  than  those  measured  by 
Hadley  (1976). 


One  way  that  we  can  generate  the  observed  nonlinearity  in  both  the  failure  envelope  and  stress-strain  curves  is  to 
account  for  the  multimineralic  composition  of  Westerly  granite.  If  the  yield  stress  in  the  Feldspars  and  micas  is 
significantly  less  than  the  yield  stress  of  quartz,  then  nonlinear  flow  of  the  weak  minerals  will  transfer  stress  to  the 
quartz,  which  will  carry  a  stress  well  above  the  average  macroscopic  loading  stress.  The  quartz  grains  fail  at  the 
higher  stress  levels  predicted  by  brittle  damage  mechanics,  but  the  apparent  macroscopic  failure  stress,  being  an 
average  of  the  flow-stress  of  the  weaker  minerals  and  the  brittle  strength  of  quartz,  is  much  smaller,  as  observed. 


A  simple  model  that  can  produce  the  required  curvature  in  both  the  failure  envelope  and  stress-strain  curve  is  shown 
in  Figure  7.  Although  we  have  not  yet  fully  explored  this  model  the  preliminary  results  are  encouraging.  As 
illustrated  in  Figure  8,  the  observ  ed  failure  stress  lies  between  the  flow  stress  of  the  weaker  minerals  and  the  elastic 
failure  of  quartz.  Moreover,  the  extra  elastic  strain  in  the  stronger  mineral  as  it  is  loaded  by  the  flow  of  the  weaker 
ones  is  close  to  the  observed  nonlinear  stain  at  failure  in  Figure  4. 
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Figure  7.  Simple  model  of  a  two-mineral  rock 
loaded  in  compression  at  a  constant 
loading  rate  e  =  constant .  The  loading 
stress  a  is  the  average  of  the  stresses 
carried  by  each  mineral  weighted 
according  to  its  eross-seetional  area.  For 
the  stronger  mineral  1, 
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with  the  same  expressions  for  the  elastic  and 
plastic  strains  in  the  weaker  mineral  2. 
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Figure  8.  Failure  envelope  for  Westerly 

granite.  The  open  circles  are  the 
strength  predicted  by  the  brittle 
damage  mechanics.  Note  that 
brittle  failure  is  an  adequate 
description  to  a  confining  stress  of 
about  -100  MPa.  The  solid  curve  is 
the  failure  envelope  measured  by 
Lockncr  (1998).  The  dashed  arrow 
shows  a  confining  stress  of  -200 
MPa.  Note  that  the  measued  failure 
stress  is  an  average  of  the  brittle 
failure  stress  predicted  by  damage 
mechanics  and  the  flow  stress  of  the 
weaker  minerals. 


The  incorporation  of  damage  mechanics  into  the  ABAQUS  dynamic  finite  element  code:  the  numerical 
simulation  of  an  underground  explosion  and  the  generation  of  S  waves  from  the  damage 


Numerical  calculation  of  this  problem  was  done  using  a  commercial  Finite  Element  Package,  ABAQUS.  A 
user-defined  subroutine,  VUMAT,  was  written  that  incorporates  the  micro-mechanics  based  damage  constitutive 
law  described  above.  The  model  dimensions  are  scaled  by  the  radius  of  the  explosion  hole,  a.  The  size  of  the  domain 
of  computation  is  / 20a  x  / 20a  (see  Figure  9).  This  domain  was  finely  meshed  with  about  120000  triangular 
elements.  The  boundary  conditions  applied  were  a  constant  uniaxial  or  biaxial  loading.  The  initial  condition  for  the 
problem  was  a  sudden  release  of  pressure  along  the  inner  periphery  of  the  hole  simulating  an  explosion.  The  entire 
calculation  was  done  under  the  assumption  of  plane-stress  conditions. 


Figure  9  Geometry 
and  mesh  design 
for  the  2D 
ABAQUS 
simulation  of  and 
explosion  shown  in 
Figure  10. 
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Figure  10.  Velocity  vectors 
and  magnitude  showing 
the  development  of  S 
waves  in  the  damaged 
area,  which  then 
propagate  to  the  far 
field.  These  calculations 
used  the  previous  D&E 
formulation  where  all 
initial  flaws  are  the  same 
size  and  are  parallel. 


The  experimental  simulation  of  an  undergroud  explosion  using  hypervelocity  impact 


Our  experiments  reported  last  year  used  a  wire  that  was  exploded  by  a  high-voltage  pulse.  When  this  explosion 
occurred  in  the  center  of  a  damaged  region  under  non-hydrostatic  load,  strong  S  waves  were  observed  to  emerge. 
However,  the  explosion  was  not  sufficiently  strong  to  nucleate  new  damage  from  the  existing  cracks.  One  was  to 
produce  a  more  intense  source  (in  2D  experiments)  is  by  using  a  high-velocity  impactor. 

Figure  1 1  shows  the  pattern  of  radial  and  circumferential  fractures  generated  in  a  Homalite  plate  by  a  high-velocity 
impactor.  Note  the  similarity  to  fracture  patterns  produced  by  underground  explosions.  We  are  currently  developing 
the  experimental  capability  to  perform  this  experiment  on  loaded  plates  and  to  measure  the  resultant  S-wave 
radiation. 


Figure  1 1.  Fracture  pattern 
generated  by  a  high-veloeity 
impact  on  an  undamaged 
Homalite  plate.  Note  that  the 
pattern  of  radial  and 
circumferential  fractures  is 
the  same  as  that  observed  in 
underground  explosions. 

The  impactor  was  an 
equidimensional  nylon  slug, 
(L/D~  1,D=0.070,;)  with  an 
impact  velocity  of  6.48  km/s 
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CONCLUSIONS  AND  RECOMMENDATIONS 

Our  primary  achievement  this  year  has  been  the  following  significant  improvements  in  the  Ashby  and  Sammis 
(1990)  damage  meehanies: 

1)  Wc  have  derived  an  expression  for  changes  in  the  internal  energy  density  associated  with  the  growth  of 
damage.  Although  we  followed  the  method  developed  by  Deshpande  and  Evans  (2008),  we  eorreeted 
fundamental  errors  in  their  formulation. 

2)  We  extended  the  previous  formulations  to  allow  non-parallel  initial  fractures. 

3)  We  extended  the  previous  formulations  to  allow  a  range  of  initial  flaw  sizes. 

Even  with  these  improvements,  we  were  unable  to  reproduce  the  observed  failure  envelope  and  stress  strain  eurves 
for  Westerly  granite  for  eonfining  pressures  above  about  100  MPa.  We  eonelude  that  a  simulation  of  the 
deformation  of  granite  at  high  eonfining  stress  requires  that  we  allow  for  plastic  flow  of  the  weaker  constituent 
minerals. 

We  have  built  our  damage  mechanics  into  the  ABAQUS  finite  element  code  and  simulated  a  2D  explosion  in 
damaged  Homalite  that  generated  significant  S  wave  energy  in  the  damage  material.  We  are  currently  ineluding  the 
above  improvements  in  the  code,  after  which  we  will  continue  with  the  simulations  of  explosions. 

Finally,  we  have  begun  2-D  experiments  that  use  a  high-veloeity  impaet  to  simulate  an  explosion.  The  higher  energy 
density  in  the  impaet  does  far  more  damage  than  does  the  exploding  wire  used  in  our  prior  experiments,  whieh 
should  allow  us  to  better  explore  the  seismie  radiation  generated  by  the  growth  of  new  damage. 
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